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Abstract 



Most hydrodynamical calculations used in heavy-ion physics ignore the effect 
of freeze-out matter carrying energy and momentum away from the expand- 
ing fluid. In a simple one-dimensional model we compare calculated energy 
density and velocity profiles, with and without interaction between fluid-like 
and freeze-out parts of the system, in order to estimate the importance of this 
effect. 
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I. INTRODUCTION 



With new data from CERN's Superconducting Proton Synchrotron (SPS) and with the 
expected completion of the Relativistic Heavy Ion CoUider (RHIC) in Brookhaven in 1999, 
nuclear physicists will be able to look at higher temperatures in larger space-time volumes 
than ever before. Of particular interest is the identification of signatures of the transi- 
tion to the high-energy phase of nuclear matter, commonly referred to as the quark-gluon 
plasma (QGP). It should be easier to measure the properties of such a transition with heav- 
ier ions because of the larger times and volumes occupied by the resulting hot and dense 
matter. With larger amounts of matter produced, we can hope that a description in terms 
of smoothly- varying locally equilibrated quantities such as temperature and fluid-dynamical 
velocity is reasonable, making hydrodynamical models especially valuable in this scenario. 
The point of view where the physics is described in terms of such average properties rather 
than in terms of individual particle motions has obvious advantages. The savings in com- 
puter simulation time allow research into a larger parameter space and a wider variety of 
initial conditions can be explored. Particle production calculations are routinely carried 
out in the framework of such models. However, these models make a drastic approxima- 
tion when assuming a sudden change from fluid-dynamical behavior to free-streaming (the 
usual effective definition of "freeze-out"). It should be kept in mind that, since our particle 
detectors register free-streaming particles far from the interaction zone, a transition in the 
description from continuous matter variables to free particles is unavoidable if the model is 
to have any predictive power or relevance to measured data. The underlying question of 
better interfacing the hydrodynamical and kinetic descriptions at freeze-out is studied by 
several groups 

In traditional hydrodynamical models, all the matter in the system of interest is assumed 
to be a thermalized fluid which expands according to idealized hydrodynamics. Within the 
fluid lies an isothermal surface which is at the "freeze-out temperature", T/q. While the 
freeze-out surface plays no role in influencing the hydrodynamic evolution, one assumes, 
for purposes of calculating particle spectra, that matter outside the freeze-out surface is no 
longer in thermal equilibrium, i.e. particle interactions cease, and particles that cross the 
freeze-out surface are free to move unimpeded to the detectors. If, for example, lepton pairs 
are measured, the decay of the freeze-out particles can contribute a significant portion of 
the observed leptons, and their spectra have interesting properties which are sensitive to the 
hydrodynamic evolution [^. 

The trajectory of the freeze-out surface is typically found by assuming that all the fluid 
is thermalized, even at temperatures below Tfo- On the other hand, calculations of particle 
spectra assume that particles near the freeze-out surface leave the fluid. The energy and 
momentum carried away by these freeze-out particles can be sizable, and should properly be 
taken into account in a self-consistent calculation. The fact that some freeze-out particles 
may return to the fluid through the moving freeze-out surface should also be accounted 
for. There has been some recent analytical investigation of how to do this by Bugaev , 
who derived equations which describe the trajectory of the freeze-out surface. Here we 
take the point of view of evolving a system from a given point in time within the context 
of a computer simulation which includes the feedback from this dynamical freeze out. A 
Godunov-type calculation Q appears well-suited for this endeavor. In this paper we show 
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how to extend the Godunov method to handle the evolution of both thermalized fluid and 
freeze-out matter, and we compare to calculations of expanding matter assuming all ther- 
malized fluid (traditional hydrodynamics). We demonstrate the effect on a one-dimensional 
example using the simplest possible relativistic equation of state (EOS) for the thermalized 
fluid, that of massless pions. These simplifying assumptions can be relaxed in the future. 
Note that at the highest energies presently in the focus of attention, where a very large num- 
ber of pions is created, even this most simplifled study may be relevant to the interpretation 
of data. 

The paper is organized as follows. In Section ^ we describe the essence of Godunov- 
type methods as typically used today and prepare the way for the inclusion of dynamical 
freeze-out. The details of how this is done are given in Section |ITH with the technicalities 
relegated to the Appendix. Section |V| summarizes the generalized algorithm implemented 
by our computer code. Section contains the results, and we conclude in Section ^ 



II. GODUNOV METHOD 

One successful approach to solving hydrodynamics on a computer is to use a so-cafled 
"Godunov- type" method 0. In such a method, the space-time continuum is discretized 
into cells, i.e. the algorithm calculates the changes in the average properties of spatial cells 
as it makes flnite steps in time. Because conservation laws tell us how quantities (such 
as energy, momentum, and particle densities) in a cell change due to the flows of those 
quantities through the cell "walls" during a timestep, it only remains to determine what the 
flows should be to give the correct hydrodynamical behavior. 

There are two common approaches within the framework of Godunov-type methods to 
calculating the flows between thermalized cells, an attempt to exactly calculate the intercell 
flow shape, and the approach of approximating the flows by uniform regions between cells. 
Both methods assume that the cells can be treated as semi-inflnite and uniform at the 
beginning of a given timestep. This is reasonable insofar as the timestep is kept small 
enough that the developing flow patterns do not extend past the middle of the cell, i.e. 
St < 6x/2, where St and Sx are, respectively, the timestep size and cell width^. Once this is 
assumed, it follows that neither a length nor a time scale exists for the development of a flow 
pattern between two cells. Therefore a similarity pattern, i.e. a pattern that depends only 
on the dimensionless variable x = x/t, develops. This is also a reasonable approximation in 
more complex geometries as long as the similarity pattern is much smaller than the size of 
the system as a whole. 

In the most ambitious approach, one attempts to flnd the actual shape of the flow 
pattern, e.g. the energy density e(x) and the fluid velocity v{x) across the cell boundary. 
These "Riemann solutions" can be derived from the continuity equations for the relevant 
quantities (in this case, energy and momentum) and the equation of state (EOS). In the 
case of a simple EOS, analytical solutions are possible For a realistic EOS, the shape 

of the flow pattern is found by numerical integration. The alternative approach proposed 
by Harten, Lax, and van Leer, and completed by Einfeldt 0, is commonly known as the 
"HLLE" method []TU| . Here the flows between thermalized cells are approximated by uniform 
regions of length proportional to the signal velocities (sound speed relativistically added to 
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fluid velocity). HLLE saves considerable CPU time because the task of calculating the actual 
flow patterns by numerical integration is replaced by choosing approximate signal velocities, 
which by conservation laws determine the conserved quantities between the two receding 
discontinuities. Discussion of the options for and choice of the signal velocities can be found 
in Ref. [0. In the present work we follow the HLLE algorithm as described there. We 
emphasize that in any Godunov-type method the calculation of the flows between cells is at 
the heart of the fluid-dynamical problem. 

For our purposes, a Godunov-type method is especially well suited, because the Godunov 
equations (see eq. (|^) below) used to evolve a cell in time are derived directly from conser- 
vation laws. Thus the validity of these equations is not restricted to ideal fluids, but they 
are true in general. Energy-momentum conservation is expressed in terms of the continuity 
equation for the energy-momentum tensor. 



. 



(1) 



The Godunov equations are obtained by averaging the time derivatives over a cell, and for 
a system that varies in only one spatial dimension, e.g. x, take the form 
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where Sx refers to the spatial extent of the cell, and the subscripts — and + refer respectively 
to the lower and upper spatial bounds of the cell. The variable t refers to the time as 
measured in the frame of the cells, which can be considered to be fixed in space (the "cell 
frame"), with matter moving between them in the course of a simulation. From inspection 
of eqs. (|1|) and (0) it is clear that T^^ represents the flow of and T^^ is the flow of T°^. 
The energy-momentum tensor T^" is defined, in the rest frame of thermalized fluid, to be 
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i.e. is the energy density e, is the momentum density, and T^^ is the pressure p 
(flow of momentum). When the fluid is given a boost v in the x-direction, we have 



E = T"" = (e + p)Y - P , 
M = = (e + p)v-f'^ , 
P = T^^ = (e + py-f^ + p 



(4) 



where 7^ = 1/(1— t;^), and we have given mnemonic names (reflecting the physical quantities 
to which these tensor elements reduce in the rest frame) E (energy density), M (momentum 
density), and P (pressure or flow of momentum density M) to the three tensor elements 

m 



of interest. {E and M are used in Ref. 



P is introduced in this work.) For concise 
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treatment, we introduce the notation Ui : E,M,P for i = 1,2,3. We also use qi : e,p 
{i = 1,2) for the energy density and pressure in the rest frame of thermahzed fluid. The 
notation F{Ui) for the flow of quantity will also be used occasionally. 

Recall that in a non-curved geometry the flows are similarity patterns, and thus the 
values of M{x — 0) and P{x — 0) remain constant throughout the timestep. It then follows 
that 



where (T^^) signifies the average value of T^'^ in the cell. This is the form that the computer 
code actually uses to evolve a cell from one timestep to the next. In the examples presented in 
this work we use a simple ultra-relativistic equation of state, so particle number conservation 
is ignored and energy- momentum conservation as presented above is sufficient. 

III. GENERALIZED GODUNOV METHOD INCLUDING FREEZE-OUT 

We wish to extend the Godunov method to handle freeze-out as well as thermahzed fluid. 
Strictly speaking, we are not extending the Godunov method per se, but simply using it in 
a way that exploits its generality. The general problem as we have mapped it out consists 
of solving for flows in three cases: between two cells of thermahzed fluid, between two cells 
of non-thermalized (freeze-out) matter, and between one cell of each kind. For this last 
case, the freeze-out surface lies somewhere in one of the two cells, and for simplicity, we 
approximate the flow with the flow due to thermahzed matter (approximation with non- 
thermalized matter is an equally valid choice). This can be done because we expect the 
density gradient to be well-behaved in the region connecting thermahzed fluid with freeze- 
out matter, as justifled a posteriori by our results (Sect. 5). 

The physics of freeze-out can be described in terms of matter that was originally ther- 
mahzed and treated as a fluid expanding to such a low density that large mean-free paths 
make fluid behavior an invalid approximation. In fluid-dynamical models of nuclear colli- 
sions, such as ours, this is assumed to happen at a well-deflned "freeze-out temperature" 
below which the particles which comprise the matter lose thermal contact. In relativistic 
nuclear collisions T/-,, ~ 100 — 150 MeV, and the freeze-out particles are expected to reach 
the detectors without any further interaction. In freeze-out matter it no longer makes sense 
to calculate flow patterns based on fluid behavior, or even to assume that local quantities 
such as energy density (e) and pressure (p) are related by an equation of state. However, 
since the Godunov equations do nothing more than express the conservation of energy and 
momentum in a discretized system, they can still be applied, as long as the correct values 
for the flows, F{E) = M and F{M) = P, (and in our case, F{P), the flow of P) are used. 
It should be pointed out that there is no relativistically covariant continuity equation for 
P, meaning that pressure is not conserved like energy and momentum. However, in weakly 
interacting matter, such as freeze-out, one can still calculate the contribution to the pressure 




(5) 



A(M) 




(6) 



5 



due to a particle of given momentum, so one can calculate the "flow of P" by considering the 
rate at which particles cross the cell wall. Thus, in the case of freeze-out matter, we assume 
A (P) = {6t/6x) [F{P)_ - F(P)+] (the calculation of F{P) is described in the Appendix), 
and for thermalized fluid we calculate P from E, M, and the EOS. 

In order to calculate flows of freeze-out matter, it is necessary to keep track of more than 
E and M. This is because for non-thermalized matter, which does not follow an equation 
of state, P does not have an implicit relationship to E and M. In principle, one should 
keep track of a number of degrees of freedom of the order of the number of particles in the 
matter. However, including just tensor element P (in addition to the usual E and M) in 
the description captures the essence of freeze-out behavior, which is that the free- streaming 
particles tend to move at a single velocity, making Pq (the value of P in the zero-momentum 
frame, where M=0, referred to as the "matter rest frame") smaller than would be expected 
for thermalized matter with the same Eq = e. Thus our calculation can be thought of as a 
first order approximation of the coupling effect. 

If the potential energy is negligible compared to the kinetic and rest mass energies, e.g. 
in a non-interacting pion gas, it is natural to express the elements of the energy-momentum 
tensor in terms of the distributions of particles [pJ| : 



T'^ = I I ^ ) , (7) 




where / (p^/T) describes the distribution of the particles. For thermalized particles viewed 
in the matter rest frame, as defined above, the function f{p^/T) is the Bose-Einstein or 
Fermi-Dirac distribution (Bose-Einstein for the pion gas considered here). The momentum 
distribution of particles immediately inside the freeze-out surface is known, because we have 
made the assumption of thermal equilibrium everywhere inside the freeze-out surface. Be- 
cause we assume a sharp freeze-out surface, the particles crossing the surface will retain the 
momenta they had while still in the thermalized fluid. We will use (|^ for the calculation 
of the elements of the energy-momentum tensor in non-thermalized matter. The only re- 
maining freedom is in specifying the distribution function, which is described in detail in 
the Appendix. 

We assume massless pions for the freeze-out particles. This is because pions are the 
lightest, and therefore most abundantly produced hadrons, so our model is closely linked to 
actual experimental data at the highest beam energies available. Massive particles could be 
used at a price of some complication of the EOS for thermalized fluid and of the calculation 
of velocity distributions. Here we limit ourselves to an extremely simple case to illustrate the 
method, which we expect to find more general applications (e.g. the nonrelativistic scenario 
of evaporation from a fluid). 



IV. STRUCTURE OF THE ALGORITHM 

Here we summarize the tasks performed by the generalized Godunov algorithm for the 
above simplified scenario at each timestep. Before evolving the system in time, the cells 
are first initialized with values for v, the cell frame matter velocities, and qi, the rest frame 
thermodynamic quantities (e and p). After all the cells are initialized, the cell frame tensor 
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elements Ui {E, M and P) are calculated. The piecewise-constant distribution is replaced by 
a piecewise-linear distribution, the slopes being calculated from the values of Ui in neighbor- 
ing cells. Thus two distinct values of f/, that correspond to two adjacent cells are associated 
with each boundary (see for details). The difference between these two Ui determines a 
flow at a given boundary. These flows give rise to second order accuracy in time by allowing 
the calculation of half-timestep quantities, including the signal velocities and flows used in 
the HLLE scheme. The new second-order flows determine the evolution of Ui inside each 
cell (throughout the code, if the Ui are changed, the corresponding and v are updated 
immediately afterward). Finally, e and v for each cell are output. At this point one time 
step is completed. 

For cells with freeze-out matter, we depart from the standard HLLE scheme, as described 
1^. If a cell is determined to be below the freeze-out temperature (based on the rest- 
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frame energy density e), the flow is calculated assuming free-streaming particles. To do this, 
we first need to parametrize the velocity distribution in terms of the tensor elements Ui, as 
described in the Appendix. To make a calculation with the freeze-out coupling turned off, 
one simply uses the standard HLLE to calculate all the flows regardless of energy density. 



V. RESULTS 

We ran most of our simulations with 6t=0.03 fm for the timestep and 6x=0.1 fm for 
the cell size. The dependence on the timestep was tested, and we found that decreasing 
6t to 0.01 fm did not alter the results appreciably. We use the equation of state of a 
thermalized massless boson gas (pions). The initial temperature was Tq = 200 MeV for 
-2.7 fm < x < 2.7 fm, and vacuum everywhere else. Initial velocities were all set to zero. 
Because of the symmetry about x = 0, only one side of the system {x > 0) is shown in the 
figures. The freeze-out temperature was chosen to be Tfo = 140 MeV, close to the mass 
of a pion. (This appears to be a realistic estimate of the freeze-out temperature p.) The 
corresponding isotherm is called the freeze-out isotherm and defines the freeze-out surface. 
In a calculation of particle production, particles that cross the freeze-out surface are assumed 
to lose thermal contact and move freely to the detectors. Therefore, the trajectory of the 
freeze-out surface is important, and is our primary concern in this work. The results shown 
were obtained by a second order calculation in time, as described in Section 4. For our 
simple problem, with a "well-behaved" EOS, we get nearly identical results for first and 
second order simulations; for equations of state with a phase transition, or other systems 
where macroscopic discontinuities may develop, a second order treatment is necessary to get 
accurate results. 

Figures and |] are contour plots of e in x-t space for the evolution of a finite one- 
dimensional slab expanding into the vacuum in both directions; only one side, x > 0, is 
shown, without and with feedback from freeze-out, respectively. The isotherms, starting 
from the right side of the diagram, represent energy densities from 5 to 70 MeV/fm^, in 
increments of 5 MeV/fm^. The thick line is the e=49 MeV/fm^ isotherm, which corresponds 
to the freeze-out temperature {Tfo = 140 MeV), or freeze-out surface. We see that it has 
a distinctly different slope in x-t space when we include the coupling between thermalized 
fluid and freeze-out. The increased inward velocity of the freeze-out surface would result 
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in more freeze-out particles being measured, because the number is sensitive to the relative 
velocity between the matter and the freeze-out surface. 

The contour plot for e also demonstrates that there are large regions of small density 
gradient in the freeze-out matter in the case of freeze-out coupling. Nearly the entire region 
outside the freeze-out surface is at around half the freeze-out energy density. This is in 
contrast to the case without freeze-out coupling, which displays a much smoother change of 
the energy- density profile. In the case with feedback from freeze-out, there are also much 
larger fluctuations in energy density, notably a region of lower density (a "bubble") that 
seems to have formed in the freeze-out matter around t = 5.0 fm, but disappears at around 
t = 7.0 fm. Whether this is indicative of new physics or merely a result of the uneven flow 
caused by the crudeness of our choice of velocity distributions, or the limited number (three) 
of moments of the momentum distribution we use, needs further investigation. 

Figures ^ and ^, the contour plots for v in x-t space, without and with feedback from 
freeze-out, respectively, corroborate our expectations and the results displayed in the e-plots. 
In the case of freeze-out coupling (Figure H), there are large regions outside the freeze-out 
surface that have little variation in velocity, indicating that the particles "free-stream" with 
no resistance. The contours are, from x = outward, v = 0.01, 0.1, 0.2, ■■■ 0.8, 0.9. The 
group of lines in the lower right corner is an artifact of the contour plotting due to the fact 
that fast rarefied matter abuts the vacuum, whose velocity is taken to be zero. Most of the 
freeze-out matter moves outward at a speed around 0.65c in the case when the feedback 
from the freeze-out is taken into account. The inhomogeneities that appear seem to be due 
to disturbances that originate at or near the freeze-out surface, where particles with a range 
of velocities cross into the freeze-out region. 
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FIGURES 





— — 7^" 














\V\\ 




\ 







1 2 3 4 5 

X (fm) 



FIG. 1. Lines of constant energy density, for the case without coupling between the freeze-out 
and thermalized fluid. The thick hne (e=49 MeV/fm^) is the freeze-out isotherm. 
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FIG. 2. Lines of constant energy density, for the case with couphng between the freeze-out and 
thermahzed fluid. The thick hne (e=49 MeV/fm^) is the freeze-out isotherm. 
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FIG. 3. Lines of constant velocity v, for the case without couphng between the freeze-out and 
thermahzed fluid. The thick line is the freeze-out surface. 
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FIG. 4. Lines of constant velocity v, for the case with coupling between the freeze-out and 
thermalized fluid. The thick line is the freeze-out surface. 

VI. CONCLUSIONS 

We have generalized the conventional Godunov algorithm to consistently include freeze- 
out matter in fluid-dynamical simulations. In particular, the energy and momentum carried 
away by the freeze-out particles and the feedback from those particles is accounted for in 
our generalized treatment. In the present paper we illustrated this dynamical freeze-out in a 
simple one-dimensional example. The net result of the dynamical treatment is a rapid eating 
into the thermalized fluid as shown by the trajectory of the freeze-out surface, as compared 
to the treatment without dynamical freeze-out. Also, there is a more sudden change across 
the freeze-out surface to freeze-out matter that is moving more or less uniformly, as expected 
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from free-streaming particles. In addition, the treatment with freeze-out couphng displays 
more fluctuations, partly due to our simphfying assumptions. 

For a more reahstic simulation we intend to use the generalized Godunov algorithm in 
more complicated geometries and with more complicated equations of state. We expect the 
method to flnd many applications. One of these may be to combine the improved freeze- 
out modeling with three-fluid hydrodynamics |jl2|. Such a combination would be especially 
useful for calculating production of particles that are sensitive to freeze-out (e.g. dileptons). 
One suspects that the effect of the freeze-out/fluid coupling on the trajectory of the freeze- 
out hypersurface will influence pion rapidity distributions at least as much as differences 
between isothermal and isochronous freeze-out [|^]. Photon spectra are also expected to be 
modifled, albeit to a lesser extent. 

We believe that fluid-dynamical calculations with an improved treatment of freeze-out 
will be very competitive with transport-theoretical modeling. This is partly because of 
great savings in CPU time and partly due to the physics insight afforded by a treatment 
dealing with average quantities, as opposed to running computer simulations with a very 
large number of chance events. 



VII. APPENDIX 

In our approach each cell is characterized by the three tensor elements E, M and P. 
These are Lorentz-transformed into the frame where M, the momentum density, is zero. 
This is referred to as the zero-momentum frame, or "matter rest frame", and we will denote 
the tensor elements in this frame by Eq, Mq = and Pq. The velocity of this frame relative 
to the frame of the cells is the "matter velocity" v. In the matter rest frame, the cells are 
characterized by Eq and Pq. For weakly interacting matter in general, Eo = e, the energy 
density, and Pq = p, the pressure. For the simple (ultra-relativistic) equation of state of our 
thermalized fluid, Pq = Eo/3. For freeze-out matter there is no relationship between Eq and 
Pq except for constraints imposed by causality. 

For thermalized massless particles, the momentum distribution is isotropic, and = p'^^, 
where we use the notation — 1 < ^ < 1 for the dimensionless velocity (in terms of the speed 
of light) along the x-axis, and p° is the time-like component of the 4-momentum. Therefore 
we can factor each integral in (^ into an integral over total energy and an integral over 
velocity in the x-direction, e.g. 



En 



(8) 

where h{^), the energy distribution with respect to velocity, is obviously a constant Eq/2 due 
to the fact that the integral over p^ is independent of C, for thermalized massless particles. 
(Note that here the subscript on Eq refers to the matter frame, while the superscripts on 
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p refer to components of the 4- momentum.) Due to this factorization for massless particles 
we can recast the momentum distribution as a velocity distribution along the x-axis. For 
other equations of state, h{^) would be defined, but not necessarily a constant. In a similar 
manner, we obtain: 

Mo = |^'rfeeMO = o, (9) 

For non-thermalized matter, we parametrize h{^) subject to the conditions that it be (1) 
even with respect to ^ = 0, which guarantees Mq = 0, and (2) non-negative. For cells near 
thermal equilibrium, i.e. Pq ~ Eq/3, we want the simplest function that differs gently from 
a constant and still satisfies the above criteria; the choice h{^) = a^^ + b {a and b constant) 
works well, where a = (45Po - l5Eo)/8 and b = {9Eo - 15Po)/8. 

For much of the freeze-out matter, Pq is too small (namely, Pq < Eg = Eq/5) to satisfy 
h{0 — + b > for all — 1 < < +1. In this case we interpolate between a function of 
the form a^^ + b and a delta function, {Eq/2)6{^). The use of the delta function reflects our 
assumption that the particles in freeze-out matter tend to move along together at nearly the 
same velocity. Occasionally it is necessary, for large Pq (> El = 3Eo/5), to use a double 
delta function of the form /i(^) = {EQ/2)[6{—d) + 6{d)], where d is a velocity. In this case 
the matter is described better by a superposition of movement at a couple velocities than 
by homogeneous particle movement. 

To calculate the flows of matter over the right-hand (left-hand) cell boundary, we calcu- 
late the integrals (^ and (|^) for particles with cell-frame velocities v{C,) greater than (less 
than) zero, where v{^) = (^ + f )/(! + ). For example, the flow of energy F{E) (due only 
to matter in the cell being considered) over the right-hand side of the cell is 

F{E)=M[ = d^^ (10) 

= Jl' d^j'-^ (vEo + {v' + l)Mo + vPo) 

where is the momentum density, viewed in the cell frame, due to particles with velocity 
between and 1. The derivatives of the rest frame quantities are simply 

dEo 



d^ 
dMo 

~df 
dPo 



HO, (11) 



di 

The evaluation of the flow of P requires the numerical integration of 



F{P)= £'^e^(0^ (12) 



d^ 

di v{i) {v^E^ + 2vM^ + Po) 
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where 7^ = 1/(1 — t'^) is the gamma factor due to the relative velocity v between the rest 
frame and the cell frame. 

Flows from two neighboring cells are added together to determine the total flow across 
the shared interface. These flows are then used by the code to update Ui at the next full 
timestep. 
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NOTES 

1. We use the convention h = c = 1. 
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